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ABSTRACT 

We develop a method to extract the "effective equation of state" of the intergalactic 
medium from the Doppler (b) parameter distribution of the low-density Lya forest. We 
C^ ■ test the method on numerical simulations and then apply it to published observations 

-Y-s ' of the Lya forest at redshifts z ~ to 4. We find that the effective equation of state 

is close to isothermal at redshift z ~ 3, indicating that a second reheating of the IGM 
took place at z ~ 3. This reheating can plausibly be identified with the reionization of 
O ! He 11 observed to occur at z ~ 3. 

Q_|. Subject headings: cosmology: theory — intergalactic medium — equation of state 
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Introduction 



It is widely believed that density and temperature fluctuations in the intergalactic medium 
5^ I (IGM) are responsible for the Lya forest of absorption lines observed in the spectra of QSOs. 

This hypothesis is supported by the results of numerical simulations that are able to reproduce 
accurately the statistical properties and main features of the observed high-resolution spectra (cf., 
Weinberg et al. 1999 and references therein). At the same time, observations of the high redshift 
Lya forest with Keck telescopes (Hu et al. 1995; Lu et al. 1996; Kirkman & Tytler 1997) and of the 
local Lya forest with the Hubble Space Telescope (HST) (Shull, Penton, & Stocke 1999; Penton, 
Stocke, & Shull 1999) have sufficient spectral resolution and signal-to-noise to permit a detailed 
study of the thermal properties of the IGM and comparisons with the output of cosmological 
simulations. 

In hierarchical dark matter dominated cosmological models, the formation of the Lya forest 
depends on two loosely related pieces of information: the evolution of the dark matter and the 
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thermal history of the IGMJ^ Recent work by Croft et al. (1999) in large part solved the problem 
of recovering the power spectrum of the dark matter from the Lya absorption spectra. Thus, 
the thermal history of the IGM remains the only piece of information still unrecovered from the 
observational data on the Lya forest. This paper complements the work of Croft et al. (1999) by 
attempting to recover the thermal history of the IGM from the data. 

It has been shown (Hui & Gnedin 1997) that in sufficiently low density IGM (cosmic 
overdensity 6^5 — 10) the gas temperature is tightly (to within a few percent) related to the gas 
density, with a power-law relation of the form T = Tq{1 + 5)'^"^, which we will call an "effective 
equation of state" because it appears to relate the density and the temperature of the gas.0 The 
thermal history of the IGM can then be described with high precision by the time evolution of the 
two parameters, Tq and 7. In the most simple physical model of the photoionized IGM adopted 
in all currently existing numerical simulations, this evolution is determined by the balance of 
adiabatic cooling due to the expansion of the universe and photoionization heating of hydrogen and 
helium. However, it is also possible that Compton heating from hard X-rays (Madau &: Efstathiou 
1999), radiative transfer effects, and time-dependent Hell reionization can be important processes 
affecting the thermal history of the IGM (Abel & Haehnelt 1999; Ferrara & Giallongo 1999). 

The most readily available piece of observational data is the joint distribution of the column 
densities and Doppler (6) parameters of absorption lines. Since the column density of an absorption 
line strongly correlates with the density of the gas in which the line originates (Miralda-Escude 
et al. 1996; Hernquist et al. 1996), we can use the column density of a line as a measure of the gas 
density. At the same time, the width of the line contains information about the gas temperature. 
Thus, we can recover the information about the evolution of the "effective equation of state" of 
the IGM from the b — N distribution of the Lya forest. This was also noted independently by 
Schaye et al. (1999). 

Undoubtedly, a better way to recover the "effective equation of state" could be developed, 
using the Hi absorption spectra rather than b — N distributions. However, this method is not 
yet developed. In addition, b — N distributions are readily available as published data, whereas 
the absorption spectra are not. Thus, we proceed with analyzing b — N distributions with the 
understanding that our results could in principle be improved upon by using the full information 
available in the absorption spectra. 

At low redshifts (z ~ 0) the amount of the intergalactic gas that is shock heated is thought to 
be larger than at higher redshifts (Cen & Ostriker 1999). This, in principle, can induce unrecovable 
systematic effects into our method. However, Dave et al. (1999) have shown that the "effective 



^The lack of understanding of how star formation proceeds prevents us from making a (in principle existing) direct 
connection between the evolution of the dark matter and the evolution of the ionizing radiation in the universe, the 
latter being the primary source for the evolution of the thermal state of the IGM. The QSO contribution to the 
ionizing radiation background that dominates up to redshift z ~ 4 is, instead, better known. 

^The true equation of state is, of course, the ideal gas equation of state P — nktT. 



DRAFT: February 1, 2008 3 

equation of state" applies down to redshift z ~ even if the fraction of gas that is shock heated 
and that does not belong to the simple power-law relation increases at low redshifts. A relation 
between the column density and the overdensity of Lya clouds also holds at z = 0. We can thus 
be confident that we can measure the "effective equation of state" of the IGM at low redshifts 
as well. Even at z = 0, there still exist a population of absorbers that are not shock heated and 
which define the lower cutoff of the b — N distribution (Dave et al. 1999). The resulting Lya lines 
are distinct, narrow features in the HST spectra (Penton et al. 1999); it is this population that we 
select by our method. In other words, the approach we adopt automatically selects against the 
shock-heated gas and thus is quite suitable for measuring the "effective equation of state" even at 
low redshift. 

The paper is organized in four sections. In § |2| we describe our method to measure the 
"effective equation of state" of the IGM. We study the output of numerical simulations for several 
cosmological models in order to understand how the thermal broadening of the lines affects the 
b — N distribution. In § |3| we describe the results of the method on the synthetic spectra created 
using the simulation outputs, and in § Q we apply the method to real observations of the Lya 
forest at O^z^A. In § |5| we discuss the results and their implications on understanding the thermal 
history of the IGM. 



2. Method 

We follow three logical steps to derive the "effective equation of state". First, we analyze 
the result of several numerical simulations using the Voigt profile fitting code AUTOVP (Dave 
et al. 1997) to build synthetic b — N distributions used to check the method. We consider several 
distinct cosmological models with different values of the ionizing intensity at 1 Ryd, J21 (in units 
of 10^^^ erg cm^^ s^^ Hz^^ sr^^), and the "effective equation of state". We also analyze different 
random realizations of the same model at the same or higher resolution. We use the Hydro-PM 
approximation to model the low-density IGM, which allows us to scan a large set of cosmological 
models and random realizations with sufficient numerical resolution, a task, currently unachievable 
with the full hydrodynamic simulations. We demonstrate that the systematic errors introduced by 
the HPM simulations are completely negligible compared to other systematic and statistical errors 
present in the method and the observational data. 



We find that the thermal broadening bx = (13 km s^^)\/T/W^K of a line with peak 
overdensity 5 (we use the notation p = 'p{l + 5) where p is the baryon density and J) is the mean 
baryon density) corresponds to the value of 6 = 5m at the maximum of the b — N distribution. 
Values of b < br are possible because of the spread of the N — p distribution and because the Voigt 
profile fitting routine introduces some lines that do not correspond to density peaks. Another 
reason why the cutoff of the b — N distribution is not sharp is that errors on b and N smooth out 
the distribution. This effect is negligible for the simulated Lya lines because we pre-determine the 
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continuum flux, zero flux levels and higher signal to noise, but it is important for the observational 
data. 

Second, we find an efficient way to locate the maximum of the b — N distribution. In order 
to use all the information available, we use the Maximum Likelihood Analysis (MLA) method to 
fit a parametric two-dimensional function to the data points in the b <Si N space. This allows us 
not only to locate the maximum of the distribution but also to study its shape. We will see that 
this could be a powerful way to extract information on the physics of the low-density Lya clouds. 
We apply the MLA method to the simulated Lya lines and demonstrate that the parameters Tq 
and 7 of the "effective equation of state" lie inside the la confidence level of the values that we 
measure for all the models and random realizations. 

Third, we apply the MLA to real data sets at redshifts z = 4 — > 0, and we measure the 
location of the maximum and the shape of the distribution. We account for the errors on b and 
A^ that, as already mentioned, spread the real shape of the distribution. At this point, we use the 
results of the simulations to interpret the b values at the maximum of the distribution as due to 
thermal broadening. This gives us the relation between T and N but not the "effective equation 
of state" . We still need to know the relationship between N and p, therefore we use the results 
we obtained from the simulations in conjunction with the Density Peak Ansatz approximation 
to determine this relationship (this is our second hypothesis based on the simulation results). It 
turns out that we only need to know the column density of the absorption lines with zero peak 
overdensity, Np, to derive the "effective equation of state" within the observational accuracy. 



2.1. Simulations 

We use an approximate simulation technique called HPM, developed by Gnedin & Hui (1998), 
to model the fiuctuations in the low-density IGM (6 ^ 10). In its essence, the HPM method 
uses a simple Particle-Mesh (PM) solver modified to account for the effect of gas pressure. The 
predefined "effective equation of state" is then used to compute the gas temperature and pressure 
at every point directly from the value of the cosmic gas density at this point. Thus, there is no 
need to introduce a special equation for the gas temperature as in a full hydrodynamic solver. As 
a result, the HPM approximation is only about 25% slower (due to the overhead of computing the 
equation of state) than a simple PM solver. It is substantially faster than a full hydrodynamic 
solver (due to both fewer computations at each time step and fewer time steps), while delivering 
results that are sufficiently accurate for our purposes. 

In order to verify the accuracy of the HPM approximation, we performed Voigt profile fitting 
for the full hydrodynamic simulations and for the HPM simulation with the same initial conditions 
and resolution used in the original HPM paper (Gnedin Sz Hui 1998). Figure || shows the b — N 
distributions for the full hydrodynamic simulation and the HPM approximation plotted on top of 
each other. The two distributions agree quite well. 
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Fig. 1. — The b — N distributions for the fuU hydrodynamic simulation (crosses) and the HPM simulation 
(filled triangles) of the same SCDM cosmological model, the same initial conditions, and the same resolution. 

Figure shows the comparison between the measured 6-parameters of the individual lines 
from two simulations. While at high b values the spread is significant, mostly due to the fact 
that the Doppler profile fitting is not a uniquely defined procedure, at low b values the HPM 
and the full hydrodynamic simulation agree on a line-by-line basis to better than observational 
uncertainties. 

In order to estimate the systematic error in the position of the peak of the b — N distribution 
introduced by the HPM simulation, we have applied our method (described below) to the two 
distributions shown in Fig. |l|. We have found that the intercepts, bo, of the 6t(-/V) relations (see 
equation ^) agree within 0.5%, and the slopes, /?, agree within 3.5%. Thus, the systematic error 
in the position of the introduced by the HPM simulation is completely negligible (at least a factor 
of 10 smaller) compared to the other systematic and statistical uncertainties (as discussed below 
in more detail). 
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Fig. 2. — Comparison of b parameters for the corresponding lines measured in the HPM approximation 
versus the full hydrodynamic simulation. The low-6 end of the distribution, where the maximum is located, 
is matched better than the observational errors. 

Numerical resolution is very important for accurately modeling the Lya forest (Gnedin 1998; 
Theuns et al. 1998; Bryan et al. 1999). The HPM simulations used in this paper have a cell size of 
10h~^ kpc and resolve the filtering scale kp over which the baryonic fluctuation is smoothed, 

(Gnedin & Hui 1998) by about a factor of three. This choice of the resolution scale has also been 
confirmed as sufficient by other authors using the full hydrodynamic simulations (Theuns et al. 
1998; Bryan et al. 1999). 

We produce synthetic spectra of the Lya forest at z = 2.85 for eight of the 25 different flat 
cosmological scenarios simulated using the HPM approximation in Gnedin (1998). The first three 
models, SCDM.2A, SCDM.2D, and SCDM.2E, are three random realizations of the same model, 
SCDM.2L is the same cosmological model with a different "effective equation of state", LCDM.3A 



DRAFT: February 1, 2008 



Table 1. Cosmological Models at z = 2.85 



Model 


fio 


!^6 


fiA 


0.,, 


h 


0-8 


n 


7 


To, 4 


ap 


J21 


(a) SCDM.2A 




0.05 








0.5 




0.7 


1.49 


0.91 


2.09 


0.17 


(b) SCDM.2D'' 




0.05 








0.5 




0.7 


1.49 


0.91 


2.09 


0.17 


(c) SCDM.2E'' 




0.05 








0.5 




0.7 


1.49 


0.91 


2.09 


0.17 


(d) SCDM.L2A'' 




0.05 








0.5 




0.7 


1.49 


0.91 


2.09 


0.17 


(e) SCDM.2L 




0.05 








0.5 




0.7 


1.33 


1.21 


2.26 


0.17 


(f) LCDM.3A 


0.35 


0.03 


0.65 





0.7 




1.04 


1.51 


1.02 


2.27 


0.30 


(g) CHDM.3A 


1 


0.05 





0.2 


0.5 




0.81 


1.49 


0.91 


1.41 


0.25 


(h) OCDM.IA 


0.3 


0.04 








0.7 




0.91 


1.24 


2.46 


1.75 


0.50 



^ A diflerent random realization of the SCDM.2A model. 
^ 512=* cells and the box size of 5.12tr^ Mpc. 

is a CDM model with the cosmological constant, CHDM.3A represents a class of Cold+Hot CDM 
models, and OCDM.IA is an open CDM model. All of these simulations have 256^ cells and 
thus have a box size of 2.56/i~^ Mpc. In order to test the effect of box size, we also analyze a 
larger simulation of the SCDM.2A model. This simulation, labeled SCDM.L2A, has 512^ cells 
and box size of 5.12/i~^ Mpc. As has been shown in Gnedin (1998), the 5.12h~^ Mpc size of the 
computational box is sufficient to reduce the cosmic variance in the column density distribution 
to a few percent. 

The parameters of the models are summarized in Table |l|, where we preserve the labeling 
used in Gnedin (1998). We used the AUTOVP Voigt profile fitting routine (Dave et al. 1997) to 
measure N and b for absorption lines along a thousand lines of sight in each simulation. This gives 
us about ten thousand absorption lines per simulation. We use this b — N distribution to verify 
the ability of our method to recover the assumed "effective equation of state" . 



2.2. Equation of State and b — N Distribution 

There exists a tight correlation between the intergalactic gas density and temperature 
in the low-density regime, where shock heating is not important (Hui &: Gnedin 1997). The 
density-temperature relation is well described by a power-law "effective equation of state" , 



T = Toil + 6) 



7-1 



(1) 



where Tq is the temperature of the gas at the mean cosmic density of baryons and 7 describes 
how the temperature changes with density. It is clear that, if 7 = 1, the temperature of the IGM 
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Fig. 3. — A scatter-plot (black dots) of the column density versus density at the peak of density fluctuations 
for the model LCDM.3A at z = 2.85, calculated by using the Density Peak Ansatz approximation and 
the simulations to measure the second derivative at the peak of the density perturbations (see § |3.2| ). The 
long-dashed line corresponds to the best fit using constant weighting of the lines and the solid line using 
weighting oc N^l"^ to account for the 0.5 slope of the logarithmic column density distribution. 

is constant everywhere (the isothermal case). The evolution of both Tq and 7 depends on the 
reionization history of the universe. 

Thus, measuring the "effective equation of state" would potentially allow us to uncover the 
evolution of the ionizing background in the universe, and b — N distribution of the Lya forest offers 
us a way to measure Tq and 7 as a function of redshift simply because the Doppler b parameter is 
related (albeit indirectly due to peculiar velocities) to the gas temperature, whereas the column 
density of an absorber is known to be tightly correlated with the gas density. 

To illustrate the latter point, we plot in Fig. |^ the column density of the absorption lines as a 
function of the overdensity at the peak of the line for the LCDM.3A model at z = 2.85 calculated 
by using the simulations together with the Density Peak Ansatz approximation (see § p.2[ ). In this 
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Fig. 4. — Distribution of log 5 versus logiV for the SCDM.2D model (black dots) at z = 2.85. The contour 
plot is the N — p distribution re-mapped on the b ^ N space using the "effective equation of state" and 
equation (^. This corresponds to peculiar velocities being set to zero. The solid line is the "effective 
equation of state" rewritten in terms of the thermal width of the lines as a function of the column density 
derived using equation (^ to relate the column density and the gas density. As explained in the text, we 
use only absorption lines in a small range of the column densities. The range of the column densities used 
and the number of lines are reported in Table EL 



case, a power-law. 



N = Njj{l + 6f, 



(2) 



is also a good fit to the distribution, but the correlation between N and 5 is not tight and the 
spread of the relation depends on the cosmological model (Zhang et al. 1996; Hui, Gnedin, & 
Zhang 1997). 



If we define the thermal broadening parameter. 



hx 



2kT\2 



niH 



(13kms"^)r/, 



(3) 
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Fig. 5. — Histograms of a vertical slice at z = 2.85 and N ~ lO^'^ cm"^ of the scatter plot in Fig. |[ The 
solid line represents the ^-distribution with peculiar velocities being set to zero (i.e., only thermal broadening 
b = bx) and the long-dashed line the ^-distribution. 

where T4 = T/IO^K, we can write the equation that relates the thermal width and the column 
density of an absorption line using (H) and (0) , 






where, 



/9 



7-1 

2r ' 



bo = (13 km s ^)Tq ^ 



1/2 (No 



Nr, 



(4) 

(5) 
(6) 



where To^4 = Tq/IO^K. For each simulation, the values of Tq, 7, Np and T are known (the last two 
are measured by fitting equation (y) to the N — p distribution as shown in Fig. y and as explained 
in § |3.2|). Therefore, using equation (0) we can derive b-r- In Fig. we plot the b — N distribution 



DRAFT: February 1, 2008 11 

for 11,520 absorption lines (black dots), taken from the spectra of the SCDM.2D model with 
J21 = 0.17 at redshift z = 2.8 using AUTOVP routine. The solid line is the thermal broadening 
6t defined in equation (|3|), and the contour plot shows the N — p distribution remapped on the 
b(d) N space using the "effective equation of state" and equation (|3|) , which corresponds to peculiar 
velocities being set to zero. It is evident that the solid line does not define the lower envelope of 
the b — N distribution owing to small-number statistics and to the intrinsic spread of the N — p 
distribution. However, both the solid line and the lower envelope track one another. Indeed, 
because we decided to fit the N — p distribution instead of finding its upper cutoff (the cutoff 
toward high column densities), the thermal broadening bT will correspond to the peak of the b — N 
distribution instead of its lower cutoff. In Fig. ^ we show this result more clearly by plotting a 
slice of the distribution at the mean column density N ~ 10^^ cm~^, averaged on a bin of width 
AlogiV ~ 0.1. The long-dashed line is the 6-distribution and the solid line is the 6-distribution 
with peculiar velocities being set to zero. Armed with this finding, in the next section we will seek 
an efficient way to locate the maximum of the b — N distribution and demonstrate that, within 
the errors of the measurement, we are able to recover 6-^. 

When we apply this method to the observed Lya forest in § B, we can only measure directly 
6t as a function of the column density A'^. We do not have a relation between the column density 
and overdensity in this case. The parameters of the "effective equation of state" are related to the 
measured /3 and Tq^ = {bo/13 km s~^)^ by the equations: 

7-1 = 2/3r, (7) 

To. = Tl,{gLf. (8) 

In § ^^ we find an approximate relationship for T and Np that allows us to solve the equations 
(|7|) and (^) given some cosmological parameters and the ionizing background intensity J21. 



2.3. Maximum Likelihood Analysis 

We define a likelihood function 

n 

^ = n ^- (9) 

where pi = p{bi, Ni\model) is the probability for the line i to have a width b = b^ and the column 
density A^ = Ni and n is the number of lines in the sample. If we describe the magnitude of the 
errors by a Gaussian distribution, then the probability function p° observed is related to the true 
one by the double convolution: 



„o6. (1 - r^)--' 



pf' = ^ '- — / / p{b, N\model) exp 



Abf ANf AbiANA 1 
2o-L 2o-^,i crb,iCrN,i J 1-^ 



dN db 
(10) 
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where ab^i and aN,i are the errors on bi and Ni respectively, A6j = bi — b, ANi = Ni — N and r is 
the correlation coefficient between the errors in b and N. Usually, the errors on the column density 
and Doppler parameter tend to be anticorrelated (r < 0), especially when the line is saturated, 
since the fit tends to preserve the equivalent width of the line. We adopt a parametric model for 
p(b, N) and the free parameters are determined by maximizing the likelihood function (Efstathiou 
et al. 1988). 

A simple and accurate method of estimating errors is to determine numerically the ellipsoid 
of parameter values defined by 

In £ = In Cn^ax - ^xUM), (H) 

where Xb{^) = t\/M is the /?-point of the x^ distribution with M degrees of freedom and ta 
is the confidence level wanted. The degrees of freedom are the number of lines minus the free 
parameters of the distribution. The parametric function, equation (|l3|), that we adopt to fit the 
b — N distribution, has nine free parameters, thus the ellipsoid of the errors is nine-dimensional. 
We are only interested in determining the errors on bo and /3. The rigorous way to find these 
should be to project the hyper-ellipsoid on the b^ — (5 space but this is computationally expensive. 
We use an approximate method that consists of maximizing the likelihood function with respect 
to all the other parameters for each value of 6o and /3. 

In general, the matrix of errors is not diagonal, and the errors on the parameters are 
correlated. The value of A'^o that makes errors on b and [3 uncorrelated (i.e., the ellipse has 
principal axes parallel to the Cartesian plane and the matrix of errors is diagonal) is the average 
column density of the sample of lines used in the fit. We choose Nq in order to minimize the 
correlation of the errors. The value of A'^o changes in each MLA mainly because the range of the 
column densities of the data-set changes. 

These are the main advantages of the method: (i) the MLA method has well-defined 
asymptotic error properties {e.g., Kendall &; Stuart 1961); (ii) there is no requirement to select 
bin widths; (iii) it is possible to account for errors on both b and N] (iv) we can also study the 
shape of the distribution; (v) all the data are used; (vi) there is no need to find the distribution 
cutoff that is not well defined because of artificial lines and metal lines. The disadvantage of this 
method is that it is difficult to test whether the assumed parametric form of p{b, N) is a good fit 
to the data. In this paper, we adopt this method but we have also used a classic minimum of 
Chi-square algorithm to fit the function p{b) = p{b, N = N) averaged over small intervals of the 
column densities, selecting several bin widths of A^. This preliminary study has been useful to 
choose the form of the parametric function that is a good fit to the data. 

In the literature, the b distribution is basically fitted with two types of parametric functions: 
(i) a Gaussian with a low b cut-off, {i.e., a three- free-parameter function) (Kirkman & Tyler 1997) 
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and (ii) a one-parameter function of the form (Hui & Rutledge 1999): 



dN 
lb 



oc —^ exp 
65 



64 



(12) 



where ha = [(6^ + 6^g^ + 6j)/(Z)+(To)]^ is the 6 value at the peak of the function. Here hj is the 
baryon smoothing scale, bres is the resolution width, Dj^ is the linear growth factor (D_|_ = {z + l)~^ 
in a flat cosmology) and uq is the rms linear density fluctuation approximately at the Jeans scale. 

We find that neither of these functions provides a good fit to the data. Instead, we use 
equation (|l2|) with the exponents of 6 in the power-law tail and in the exponential cutoff as free 
parameters. We notice that, as the column density increases, the power-law tail of the distribution 
becomes more narrow and the exponential cutoff becomes more sharp. We have also tried a simple 
power-law 6~" with a low 6 cutoff, motivated by the fact that the exponential cutoff is very sharp, 
but we have found that this form of the distribution is not a good fit to the data. 

The form of the parametric function that is the best fit to the simulated Lya forest and that 
we will use in the rest of the paper is. 



p{b,N) 



■ MN"" ( J 



exp 



(13) 



where, 



6* 







(14) 
(15) 



J^min and r^max 



are the low and high column density 
*dt is the Gamma function. The 6 value at the 



Here, N is the normalization constant, 

cutoffs of the distribution and T{z) = /q°° t^~^e~ 

maximum of the function is 6a/. As previously stated, we assume 6t = 6a/. We do not use the Hui 

& Rutledge relation for h^ because the Jeans smoothing scale in this expression is needed only if 

we want to relate the dark matter density to the baryon density and temperature. We already 

include this smoothing scale in the N{p) expression that we derive in § 3^. We assume that the 



broadening of the lines due to the simulation/observation resolution is negligible in the range of 
the column densities considered. We also find that the values of the parameters $i, ^2 and 932 are 
practically zero in all the MLAs for the simulations and observations as well. The shape of the 
distribution is then determined by 6 parameters: 60 and (3 which determine the value of 6a/, <^o 
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Fig. 6. — Value of the free parameters in equation ( |13| ) as a function of the correlation coefficient r. The 
parameters are normalized to their value at r = 0. 

and (fi which are related to the power law tail, $o related to the exponential cutoff, and a the 
slope of the column density distribution. 

We do not know the correlation coefficient r of the errors on 6j and Ni. However, as is shown 
in Fig. ^, when we change r between zero and -0.95, the values of the parameters determined with 
the ML A change by at most 20%, well below the observational accuracy or the accuracy of our 
method. 



3. Testing the Method 



DRAFT: February 1, 2008 



15 



Table 2. Simulated Lya forest 



Model 


number 
of lines'^ 


log iVmin 

(cm-2) 


log Nmax 

(cm-2) 


log iVo 
(cm-2) 


feo 
(kms^^) 


/3 


V>o 


VI 


$0 


a 


(a) SCDM.2A 


11500/1716 


13.2 


14.6 


13.7 


16.8 


0.18 


3.1 


3.5 


3.6 


1.9 


(b) SCDM.2D'' 


11520/1744 


13.2 


14.6 


13.7 


18.0 


0.16 


3.9 


3.7 


3.8 


1.8 


(c) SCDM.2E'' 


10222/1204 


13.2 


14.6 


13.7 


18.5 


0.15 


4.0 


3.7 


3.7 


1.9 


(d) SCDML.2A^''" 


12131/1373 


13.3 


14.5 


13.7 


19.0 


0.16 


4.1 


3.1 


2.5 


1.8 


(e) SCDM.2L 


6980/1370 


13.0 


14.9 


13.7 


19.2 


0.11 


6.9 


3.7 


2.2 


1.7 


(f) LCDM.3A 


9233/1182 


13.2 


14.4 


13.7 


19.4 


0.14 


3.6 


3.7 


4.9 


1.7 


(g) CHDM.3A 


11812/1496 


13.2 


14.5 


13.7 


17.6 


0.17 


2.8 


3.7 


4.2 


1.9 


(h) OCDM.IA 


4289/271 


13.5 


14.9 


13.7 


26.1 


0.04 


12 


0.1 


13 


1.9 



^ A difFerent random realization of the SCDM.2A model. 

"" 512^ cells and the box size of 5.12/i~^ Mpc. 

'^ Total number of lines/number of lines in the column density range Nmtn < A' < Nmax- 

3.1. The b-N Distribution 



The number of absorption lines in each simulation is about 10^, but we are restricted to 
use lines in a range of the column densities approximately one decade wide. This reduces the 
number of lines to about 10'^. The upper limit on A^ is determined for each simulation by 
imposing a maximum overdensity at the line center, 6 ~ 10, which is the limit below which HPM 
simulations produce accurate results. The lower limit on the column density is determined by the 
completeness requirement of the distribution and by intrinsic problems associated with the Voigt 
fitting routine AUTOVP. The relative errors on b and A^ are about 10%, smaller than typical 
observational errors, since we assumed a higher signal-to-noise ratio than is usually afforded by the 
observational data (in order to test our method in more challenging conditions) and because we 
do not have uncertainties due to continuum fitting and zero fiux level. Therefore, the convolution 
with Gaussian errors does not affect the shape of the distribution in a significant way. 

In Table |2| we show the values of the parameters that give the best fit to the b — N distribution. 
We note that the value of <pi is surprisingly constant for all the models but one. Also, the lack of 
lines in the region of large A^ and b turns out to be a real feature of the distribution and not an 
apparent feature due to a decreasing number of lines. The width of the b distribution at fixed A'^ 
is proportional to 1/ip and decreases with A^ as l/{^o + ^iN). The value of $, instead, varies 
between four and six and is directly related to the width of the N — p distribution. 

In Fig. we present the result with which we are more concerned in this paper: that the 
parameters of the equation of state can be measured from the b — N distribution within the 
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Fig. 7. — Error ellipses for the models at z — 2.85, listed in Table |i|. The solid lines show la error contours 
and the dashed lines show the 2cr error contours. The crosses are the expected values of bo and (3 calculated 
with equation (||). The size of the crosses corresponds to the systematic error introduced by the HPM 
approximation. 

precision of the method. Each panel shows the best fit (black dot) and the la and 2a error 
contours for bo and /3. The triangle and the cross are the theoretical values derived using equation 
(^ for two slightly different measurements of the parameters Np and F in the equation (^) . In one 
case, we fit the N — p distribution with constant weighting for all the lines (triangles), and in the 
other (crosses) with a N^'"^ weighting in order to correct for the effect of the decreasing number 



of lines with increasing the column density. For a more detailed explanation see § 3.2. In all the 



models, the expected values of &o ^-nd (3 lie near the la contour level of the measurement, as they 
should be if the errors are estimated correctly. 
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3.2. Relation Between Column Density and Overdensity 

The Density Peak Ansatz (DPA) approximation is based on the assumption that a single Lya 
line arises from a peak in the cosmic gas density along the line of sight to a distant quasar. Then, 
given the value of the cosmic density at the peak, and the second derivative of the density along 
the line of sight (the first derivative vanishes at the peak), a value for the column density can be 
associated with the peak as follows (Gnedin & Hui 1996): 



Nui = (3.35 X 10^'^cm-^)r( 



QbffA /a5\ fl + z 
''''" I 0019 j [l^iJ 



X (1+5)2-0.7(7-1) 



;i_ 0.35(7-1)) 



d^ ln(l + 5) 



and X being measured in Mpc. Here J21, is the ionizing intensity, 

/ JyOydv /v 1 



-1/2 



(16) 



21 = ^ ^-, — X 



^(Tydv/v 10-21 erg cm-^g-i Hz -"^sr-^' 

where Jy is the radiation specific intensity and a^ is the Hi photoionization cross section. 
Comparing equations (|2|) and ([T^ we find that one can write the inverse of the square root of the 
average curvature at a density peak, i.e. the typical comoving size (expressed in Mpc) of a Lya 
cloud, as a function of the overdensity, 

where H and ^ are model dependent. We hypothesize that the average curvature of a density 
peak is a function of the non-linearity of the baryonic density field. An accurate indicator of the 
non-linearity of the gas distribution of a given cosmological model is the linear density fluctuation 
at the filtering scale, ap (Gnedin & Hui 1998). We therefore fitted H and ^ as linear function of 
ap for 25 flat cosmological models (Gnedin 1998). The result of the fit is given in Table ^. The 
parameters in equation (pi) are measured using the minimum of Chi-square analysis of the N — p 
distribution with a weighting on each line oc N^'"^ in order to account for the 0.5 slope of the 
logarithm of the column density distribution. 

Combining equations (|l7| ) and (|l^), we write the relationships, 

r = 2 - 0.7(7 - 1) + C =^ 1.68 + Q.m^ap - 0.7(7 - 1), 

2 

0.019; \j2i)\ 4 ; [1 -0.35(7 -i)]V2- 

Both r and N-p are insensitive to the value of ap if we make a conservative assumption that it is 
in the range 2 it 1 (Gnedin 1998). From now on, we use a value ap = 2. 



iV,M3,35 . .0-c„,-,T- m\ {'4) f^r °-!!-"°^^-„ . m 
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Table 3. Parameters of the best fit 



X 



S = a + bap 
^ = a + bap 



0.19 ±0.02 
-0.319 ±0.02 



-0.025 ± 0.01 
0.068 ± 0.01 



0.005 
0.013 



Finally, if we define N^ = N-pT^'"^ , we can write the equations that relate the measured 



parameters (3 and Tq^ to the "effective equation of state" parameters, 

3.634/3 



7-1 



f0,4 



1± 1.4/3' 

^t I P I 



■^0,4 



iVn 



1+1.4/3 



where. 



N^ (4.7x 10i2cm"2) /sij/i2y /o.5\ /1 + 



1 + 1.4/3 



A^o A^o \^0.019y V^2i/ V 4 y y 1 + 0.128/? 

Propagating the errors on (3 and 6o) '^s write the errors on the "effective equation of state" 
parameters, 

A(7 - 1) 1 A/3 



(19) 
(20) 

(21) 



7-1 
ATo,4 

^0,4 

If m/No = (6o/13 kms^ 



1 + 1.4/3 /3 ' 



1 + 1.4/3 



A6n 



6o 



+ 



A/3 



'ln^-1.41nV, 

, iVo 13/ 1 + 1.4/? 



1 

2\ 2 



(22) 
(23) 



'1M.4 



the errors are uncor related and ATq 4 has a minimum value. 



4. Observed Evolution of the Equation of State 

4.1. Observations 

We apply our method to the published data sets of Lya absorption lines. In this paper, we 
analyze eight data sets, but the "effective equation of state" can be measured with much better 
accuracy if we use more observations. In order to obtain good results, we need to use the data 
sets with small errors on b and A^ and with as many lines as possible inside the range of the 
column densities where the completeness of the sample is good. For the high-redshift Lya forest 
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Table 4. Observational data 



QSO 


Telescope 


FWHM 


S/N 


Voigt fitting 


number 






km s~^ 




code 


of lines 


QOOOO-26 


Keck HIRES 


6.6 


20-80 


VPFIT 


337 


APM 08279+5255 


Keck HIRES 


6 


80 




380 


Hu et al.'' 


Keck HIRES 


8 


50 




1056 


GB 1759+7539 


Keck HIRES 


7 


25 




318 


HS 1946+7658 


Keck HIRES 


7.9 


15-100 


VPFIT 


459 


B2 1225+317 


Kitt Peak 4-m 


18 


2-20 




159 


Penton et al. 


HST/GHRS 


19 


10-42 




80 



Note. — The data for the quasars at high-z are from Lu et al. 1996 (QOOOO-26), Ellison et al. 1999 
(APM 08279+5255), Outram et al. 1999 (GB1759+7539), Kirkman & Tytler 1997 (HS 1946+7658), 
Khare et al. 1997 (B2 1225+317). 

'^Collection of high-z quasars Q0014+813, Q0302-003, Q0636+680, Q0956+122 (Hu et a/., 1995). 

•"Collection of low-z quasars: AKN120, ES0141-G55, FAIRALL9, MARK279, MARK290, 
MARK421, MARK509, MARK817, MARK501, PKS2155, Q1230+0115, H1821+643, ES0141-G55 
(Penton et al. 1999). 

we use Keck-HIRES observations of the quasars QOOOO-26 (Lu et al. 1996) and APM 08279+52550 
(Ellison et al. 1999) at z ~ 4, Q0014+813, Q0302-003, Q0636+680, Q0956+122 (Hu et al. 1995), 
GB1759+7539 (Outram et al. 1999) and HS 1946+7658 (Kirkman & Tytler 1997) at z ~ 3. At 
intermediate redshift we use the data-set of Khare et al. 1997 for the quasar B2 1225+317 at 
z ~ 1.9. Finally, we use HST observations of the local Lya forest kindly provided by Penton, 
Stocke, & Shull (1999), using the Goddard High Resolution Spectrometer (GHRS) with 19 km s~^ 
resolution. 

We use a subsample of 337 Lya absorption lines of the complete Lu et al. (1996) data set 
within a redshift range 3.425 < z < 3.977. The upper cutoff is chosen because this is where the 
UV ionizing intensity from the quasar is expected to roughly equal that of the UV background. 



^The data presented herein were obtained at the W. M. Keck Observatory, which is operated as a scientific 
partnership among the Cahfornia Institute of Technoiogy, the University of California and the National Aeronautics 
and Space Administration. The Observatory was made possible by the generous financial support of the W. M. Keck 
Foundation. 
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2.5 





Fig. 8a 



Fig. 8b 



Fig. 8. — (a) Distribution of log 6 versus log TV for Lya lines for QOOOO-26 at (z) ~ 3.7, with errorbars from 
VPFIT. (b) Comparison between the b distribution for QOOOO-26 (long-dashed hne) and for model LCDM.3A 
(solid line) at z = 2.8 and N ~ 10^'^ cm^^. The large errors on b and N of the observed Lya forest spread 
considerably the true distribution. 

The data set from Hu et al. (1995) is a collection of four quasars and consists of around 10^ 
lines in the redshift range 3.425 < z < 3.977. Unfortunately, the published list of lines does not 
have errors on the h parameter and the column density. In order to use this large list of lines, 
we decided to produce synthetic errors, assuming that their amplitude and distribution follows 
the same statistics as the Lu et al. (1996) sample. We chose this sample because it has similar 
resolution and S/N. Unfortunately, the Voigt profile fitting routine is not the same. 

The local Lya data set (Penton et al. 1999) consists of 79 lines from 15 Seyfert galaxies and 
quasars detected in observations with the HST/GHRS, using both pre- and post-COSTAR optics. 
We do not use the pre-COSTAR observations because this sample appears to have systematically 
larger 6- values due to spacecraft wobble (Penton et al. 1999). This reduces the sample to 43 lines, 
which are distinct Lya absorption features at > 4(7 significance, with (6) = 35.0 it 16.6 km s^^. At 
low redshift, a typical line of sight to a QSO sA. z ^ 0.05 may have fewer than 10 absorbers. As a 
result, it is not possible to have an uniform sample of lines belonging to the same line of sight such 
as for the high-redshift Lya forest. However, the ensemble of \ow-z Lya absorbers, drawn from 15 
sight-lines, adequately characterizes the distribution of 6-values down to the 19 km s~^ spectral 
resolution of the HST/GHRS. Additional Lya lines are currently being observed with the Space 
Telescope Imaging Spectrograph (STIS) and will double the sample of lines. 
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In all the samples, we do not include identified metal lines. We only use absorption lines in 
a limited range of the column densities. The lower limit is determined by imposing the following 
two conditions: 

• the completeness of the sample must be close to one; 

• the width of the b distribution must be monotonically decreasing in the range considered. 

The first condition is usually the most stringent. In the case of QOOOO-26, the two conditions are 
comparable. The upper limit on the column density is set by imposing the following constraints: 

• overdensity of the line 6^10 (range of validity of the HPM approximation). 

• N-max 1^ lO^'''^ cm~^ in order to avoid heavily saturated lines. 

By restricting ourselves to this range of the column densities, we are more confident of the 
reliability of the results from the Voigt profile fitting algorithm. Indeed, we find that the shapes 
of the ^-distributions obtained by different authors, using different Voigt profile fitting algorithms, 
are very similar in this range of the column densities. It is however possible that fluctuations of 
the ionizing background (Haardt & Madau 1996; Fardal et al. 1998) can produce different mean 
absorption of the quasar spectra and result in a spatially fluctuating value of Nj. However, this 
effect is much smaller than the errors of the measurement. We will quantify the change of the 
IGM temperature produced by a change in the ionizing background below. 

We discard lines with errors on 6 or A^ greater than 50% in order to obtain more precise 
results. Large errors degrade the precision of the flnal result. On the other hand, the number of 
lines of the sample decreases if we consider only the lines with small errors. By performing several 
MLAs, with various error thresholds, we find that limiting the relative errors to about 50% is the 
best choice. 

In Table § we summarize the main properties of the data sets used in this paper. To give an 
idea of the magnitude of the errors and how they can affect the shape of the b distribution, we 
show in Figs. ^ and^a the data set of Lu et al. (1996), with error-bars and the b distribution at 
N ~ 10^^ cm~^ compared to a typical simulated b distribution. 



4.2. Results 

In Table || we report the parameters for the best fit from the MLA. We use the same form of 



the function as adopted in the simulations (see equation [13|). Fig. ^ shows the values of 5o and /3 
and the error contours at Icr and 2(T confidence levels. Each panel shows the results for a different 
quasar or data set. In Table ^ the quasar and the mean redshift of the Lya forest are reported 



along with the label that appears on the upper right corner of each panel. Fig. IC is analogous to 
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Table 5. Observed Ly-a forest 



QSO 



(2) number log N 
of lines'^ (cm 



2 



log No 



) (c 



bo 

(km s-^) 



•fio 



Vi 



$ 



(a) QOOOO-26 


3.7 


195 


(b) APM 08279+5255 


3.42 


186 


(c) Hu et ar 


2.85 


636 


(d) Q0014+813 


2.95 


159 


(e) Q0302-003 


2.85 


160 


(f) Q0636+680 


2.75 


193 


(g) Q0956+122 


2.85 


180 


(h) GB 1759+7539 


2.72 


13.5 


(i) HS 1946+7658 


2.7 


221 


(1) B2 1225+317 


1.9 


99 


(m) Penton et al} 


0.06 


43 



13.5 


14.15 


23.6 


0.11 


3.6 


0.1 


4.6 


1.7 


13.0 


13.7 


23.1 


0.09 


2.6 


0.3 


5.6 


1.7 


13.0 


13.7 


24.3 


0.04 


4.6 


0.01 


3.1 


1.6 


13.0 


13.7 


26.4 


0.07 


5.9 


0.1 


2.1 


1.7 


13.0 


13.7 


26.9 


0.11 


4.8 


0.2 


2.8 


1.7 


12.8 


13.7 


23.9 


0.05 


4.7 


0.06 


3.3 


1.6 


12.8 


13.7 


25.7 


0.04 


4.3 


0.02 


4.4 


1.6 


13.9 


13.7 


23.2 


0.16 


3.1 


1.5 


2.0 


2.1 


13.0 


13.6 


22.4 


0.06 


3.3 


0.05 


7.2 


1.8 


13.2 


13.7 


24.4 


0.1 


6.2 


0.5 


1.2 


1.8 


12.5 


13.48 


34.2 


0.35 


12.4 


0.7 


1.4 


1.6 



='■'' See notes on Table §. 

'^ Number of lines used in the MLA in the column density range Nmin ^ N < Nmax 



Fig. ^ but shows the value and error contours for 7 — 1 and Tq 4 derived using equation (|T9|) . Finally 
in Fig. 11 we show the values and the la error contour for 7 — 1 and Tq (computed with equation 



In each panel we have the result for two different values of the ionizing intensity J21. The 
values chosen for J21 at different redshifts are based on observational (Lu et al. 1996; Shull et al. 
1999) and theoretical (Haardt & Madau 1996; Fardal et al. 1998; Valageas & Silk 1999) results. 



We adopt as minimum values of J21: J™" (2; = 4) = 0.2, J; 



21 



3) = 0.5, Jat 



2) = 0.2 



and J> 



21 



0) = 0.01. The maximum values are, J2i^^{z) = 2 x Ji 



21 



for all z. 



The error ellipses are strongly elongated and tilted with respect to the Cartesian axis. This 



happens because N^ y^ Nq. In this case, as shown in equation (23), the error on Tq depends on the 
value of 7 — 1. It is impossible to remove this correlation of the errors by changing the value of the 
normalization parameter A'^o because this would introduce exactly the same amount of correlation 
from the errors on bo and /3 due to the limited range of the column densities of the lines. Using a 
smaller value of the minimum column density of the distribution of absorption lines would be the 
only way to reduce this correlation. This means that, in order to have N^/No = 1, we would need 
a data set that is also complete for lines arising from the density fluctuations in the underdense 
regions {6 < 0). This would probably imply the use of a more complicated form of the parametric 
function if the b distribution starts to become narrower in underdense regions as can be expected 
on theoretical grounds. 

In Fig. |l^, we show the final result of this work, the evolution of the "effective equation of 
state" of the IGM with la error-bars and for two different values of the ionizing radiation J21 , 
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Fig. 9. — Error ellipses in the [3 ®hQ space for the data sets listed in Table ||. Solid lines show la error 
contours and the dashed lines show the 2a error contours. 

chosen as explained above. The results are also listed in Table g and Table for reference. The 
point at z ~ 4 is the weighted average of the observations (a) and (b), and the point at z ~ 3 is 
the weighted average of the observations (d), (e), (f), (g), (h), (i), where the letters refer to the 
list in Table ^. 

The lines in Fig. ^ show simple photoionization models for the thermal history of the 
universe computed using the method of Hui & Gnedin (1998). In all the models, hydrogen is 
reionized suddenly at z = 7. The solid line shows the model where He ii is not reionized at all, the 
short-dashed lines show models with sudden helium reionization at z = 3, the long-dashed line 
shows the same model but with the He ii photoheating rate artificially increased by a factor of four 
following the suggestion by Abel & Haehnelt (1999). Finally, the dot-dashed line shows a model 
where the He ii reionization is more gradual, and it provides a remarkably good fit to the data. 
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Fig. 10. — Error ellipses in the 7 — 1 (g) logTg space for the data sets listed in Table ||. Solid lines show la- 
error contours and the dashed lines show the 2a error contours. 



Table 6. Time evolution of the observed "effective equation of state" . 



{z) 7-1 To/lO^ K To/lO^ K ( J21 X 2) ATo/Tq J21 

1.72 
2.21 
1.57 
0.34 



3.56 ±0.27 


0.38 ±0.18 


1.98 


2.75 ±0.27 


0.22 ±0.10 


2.52 


1.90 ±0.30 


0.32 ±0.30 


1.77 


0.06 ±0.01 


0.85 ±0.30 


0.47 



16% 


0.2 


9% 


0.5 


32% 


0.2 


37% 


0.01 
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Fig. 11. — 1(7 error ellipses in the 7 — 1 (g) log To space for the data sets listed in Table g. The two ellipses in 
each panel are calculated using J21 = >/2i*"i J^i"^- Values are reported in the text. 

4.3. Sanity Check 

We would like to emphasize that the thermal history of the universe which we derive from 
the observational data is quite different from the one assumed in all the simulations we used 
to test the method: the temperature is significantly higher and the "effective equation of state" 
flattens out at z = 3. As a sanity check, and in order to test the reliability of the method for the 
thermal history similar to the one observed, we ran a simulation of a "realistic" LCDM model 
with 0,Q = 0.3, r^A = 0.7, r^b = 0.04, h = 0.7, n = 1, cg = 0.91 and an equation of state shown by 
the long-dashed line in Fig. ^, which is sufficiently similar to the observational data. We follow 
exactly the same procedure used to analyze the observations to derive Tq and 7. 

Fig. O demonstrates that our method works in this case as well, despite the fact that it was 
developed, based on the simulations with the thermal history different from the one observed. In 
that figure, we show the "realistic model" over mentioned at redshifts z = 3.69,3.34,2.85,2.31, as 
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Fig. 12. — The evolution of 7 (a) and To (6) with redshift. The point with error bars at z ~ 4 is the weighted 
average of the observations (a) and (b), and the point at z ^ 3 is the weighted average of the observations 
(d), (e), (f), (g), (h), (i), where the letters refer to the list in Table 13. The solid line represents the theoretical 
model with a sudden Hi reionization at Zrd = 7 and with Tjoi ~ 2.5 x 10^ K. The dashed lines show the 
effect of sudden He 11 reionization in the optically thin approximation (short- dashed line) and with the Hen 
photoheating rate quadrupled (long-dashed line). The model shown by the dot-dashed line has a quadrupled 
He II heating rate and a more gradual He 11 reionization. 

well as the assumed thermal history (dashed line). One can see that our method reproduces the 
correct answer within the errors. 

Comparing Table || with Table |5|, we note that the simulated and observed b — N distributions 
are very similar apart from the values of the parameter ipi. In all simulations but one, the 
parameter ipi, which measures how much the ^-distribution narrows at high column densities, is 
about 3.7. In the observational data, it is indistinguishable from zero, which implies that the 
observed 6-distribution has the same width for all values of the column density. In order to 
understand this difference between the 6-distributions we applied our method to a sample of 60 
Lya lines at z = 3 from a cosmological hydrodynamic simulation of Miralda-Escude et al. (1996), 
where gas shock heating is included. For this simulation we found a value (^1 = 0.1, consistent with 
the observations. Our interpretation of this result is that shock heating spreads the 6-distribution 
at high column densities. This does not affect in any way the validity of our determination of the 
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Table 7. Measured "effective equation of state" of the observations at z ~ 4 and z ~ 3. 



(z) 



7 



To/lO^ K 



To/lO^ K (J21 X 2) 



ATo/To 



J21 



3.70 ± 0.27 (a) 


0.35 ±0.2 


2.04 


3.42 ±0.27 (b) 


0.47 ±0.4 


1.78 


2.95 ±0.25 (d) 


0.23 ±0.3 


2.69 


2.85 ±0.25 (e) 


0.34 ±0.3 


2.21 


2.75 ±0.25 (f) 


0.17 ±0.2 


2.43 


2.85 ±0.25 (g) 


0.14 ±0.3 


2.97 


2.72 ± 0.33 (h) 


0.30 ± 0.3 


1.77 


2.70 ±0.30 (i) 


0.14 ±0.25 


2.69 



1.78 


15% 


0.2 


1.51 


25% 


0.2 


2.47 


30% 


0.5 


1.94 


30% 


0.5 


2.27 


20% 


0.5 


2.82 


25% 


0.5 


1.58 
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Fig. 13. — The evolution 7 (a) and Tq (b) as a function of redshift for the "realistic" LCDM model with 
sudden He 11 reionization at z ~ 3 and a photoheating rate 4 times the optically thin approximation (dashed 

line). 
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"effective equation of state" , but instead can, in principle, be used to measure the fraction of the 
gas at a given density that is shock-heated. 

5. Summary and Discussion 



We note two main features of Fig. 12: the value of 7 — 1 at 2: ~ 3 is very low and the 
temperature at 2: ~ 3 is too high. These features are only marginally consistent with simple 
photoionization models in which Hi and Hell are photoionized simultaneously at z ~ 5. The 
direct interpretation of this disagreement is that H i reionization occurs quite late, after 2; = 5 or 
so. This, however, contradicts the observation of Lya emission lines from quasars and galaxies at 
z > 5 (Weymann et al. 1999; Chen et al. 1999), which would be impossible if the universe was 
neutral at z^5 (Miralda-Escude & Rees 1998). Another interpretation of our results is that the 
universe reheated a second time at 2; ~ 3. Reheating will increase the temperature of the IGM and 
flatten the effective equation of state (7 — 1 decreases) . One of the most plausible source of the 
secondary reheating is He 11 reionization. Recent observations of the He II Lya forest also support 
the conclusion that Hell reionization occurs at 2; ~ 3 (Reimers et al. 1997; Anderson et al. 1999; 
Heap et al. 1999). Abel & Haehnelt (1999) have recently shown that radiative transfer effects can 
increase the He 11 photoheating rate respect to the optically thin approximation by up to a factor 
4. This results in an increase of the temperature by a factor 1.5-2.5. 

The dashed lines in Fig. ^ show the effect of a sudden He 11 reionization at 2 = 3 in 
the optically thin approximation (the short-dashed line) and with the He 11 photoheating rate 
quadrupled (the long-dashed line). The model shown by the dot-dashed line has the quadrupled 
He II heating rate and a somewhat gradual He 11 reionization. Despite being so simple-minded, the 
latter model provides a remarkably good fit to the observed thermal history of the IGM. 

Other possible sources of heating in the IGM are photoelectric dust heating (Nath, Sethi, &: 
Shchekinov 1999) and Compton heating by hard X-rays. Recent results (Madau & Efstathiou 
1999) indicate that the second could be an important heating source at redshifts larger than 
2 ~ 10 but photoheating is dominant at lower redshifts. Photoelectric heating of dust grains could 
be a minor source of heat but some recent observations show anly a small amount of dust in the 
high-2 Lya forest (Outram et al. 1999). 

We would like to emphasize here that, while our error bars on the thermal history of the IGM 
are quite large, a substantial improvement can be obtained with a larger data set and extending 
the completeness of the sample to low-A^ lines. The error bars are still dominated by the statistical 
errors, as can be seen from comparing Fig. [^ and Fig. |^ With a larger data set, the errors on Tq 
and 7 can be reduced up to a factor of three, after which the systematic errors become significant. 
The most significant systematic error is the lack of knowledge of the intrinsic shape of the 
b — N distribution; with the larger data set this error can also be reduced by introducing a more 
complicated fit to the overall distribution. The errors on the parameters /3 and 60 measured with 
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the MLA method are proportional to the square root of the number of Unes of the sample. For a 
typical simulation with 10^ lines, we have A(3 ~ 0.06 and A6o ~ 1-5 km s~^. For the Hu et al. 
(1995) with 636 lines, we have (3 ~ 0.08 and A^g ~ 2 km s~^. The errors on the "effective equation 
of state" parameters are approximately A(7 — 1) ^ 3.634A/3 and ATq/Tq ^ 2A&o/^o- At redshift 
2; ~ 3 we have six QSOs spectra and the errors are A(7 — 1) ~ 0.1 and ATq/Tq ~ 9%. Thus, more 
observations at redshift z ~ 4 and z ~ 2 are needed to reduce the current errors. Reducing the 
errors on the Doppler parameters and the column densities of the lines are important as well, in 
order to obtain an accurate measure of the "effective equation of state" and of the shape of the 
b — N distribution. Also, more absorption lines for the local Lya forest are required in order to 
put strong constraints on the evolution of the "effective equation of state" of the IGM at low z. 

One of the major issues for future work concerns the distribution of Doppler 6- values for the 
low-z Lya absorbers (Penton et al. 1999). Although we worked here with just 43 absorption lines 
at z < 0.06, many more lines of sight are currently being analyzed with HST/STIS spectroscopic 
data. Thus, theoretical predictions of the thermal evolution of the IGM at low redshift will soon 
be testable. One of the intriguing predictions of many hydrodynamic simulations (cf., Cen & 
Ostriker 1999) is that the relative fractions of warm (photoionized) absorbers and hot (shocked) 
shift considerably between z = 2 and z = 0. The empirical distribution of b parameters for Lya 
absorbers (from HST) and Ly/3 absorbers (from FUSE) should reflect this shift and can provide a 
test of these models. 

We are grateful to Steve Penton for the as yet unpublished data set on the local Lya forest 
and to Mark Fardal for the list of published data on the high-redshift Lya forest. We thank David 
Tytler, Len Cowie, and Limin Lu for allowing us to use the published b — N distributions. We also 
thank Romeel Dave for letting us use his AUTOVP automated Voigt profile fitting software. This 
paper was significantly improved as the result of fruitful conversations with Andrew Hamilton and 
Jordi Miralda-Escude. This work was supported in part by grants from NASA (NAG5-7262) and 
NSF (AST96-17073). Calculations were performed the NCSA Origin2000 array under the grant 
AST-960015N. 
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